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Fokker-Planck simulations for rotating star clusters: 
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ABSTRACT 

We have carried out iV-body simulations for rotating star clusters with 
equal mass and compared the results with Fokker-Planck models. These 
two different approaches are found to produce fairly similar results, al- 
though there are some differences with regard to the detailed aspects. We 
confirmed the acceleration of the core collapse of a cluster due to an initial 
non-zero angular momentum and found a similar evolutionary trend in the 
central density and velocity dispersion in both simulations. The degree of 
acceleration depends on the initial angular momentum. Angular momen- 
tum is being lost from the cluster due to the evaporation of stars with a 
large angular momentum on a relaxation time scale. 
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1 INTRODUCTION 

There are two different approaches to study the dy- 
namical evolution of collisional stellar systems (for ex- 
ample, globular clusters): statistical approach and di- 
rect integration of the TV-body equations of motion. 
Among the statistical methods, the Fokker-Planck 
(henceforth referred as FP) equation, which is a sec- 
ond order approximation of the collisional Boltzmann 
equation, has been frequently used. The FP equation 
is solved by using either Monte Carlo techniques (e.g., 
the series of eight papers, from Spitzer & Hart 1971 
to Spitzer & Mathieu 1980, see also an alternative 
approach by Henon 1971, and for recent adaptations 
Giersz 1998, Joshi, Rasio, & Portegies Zwart 2000, 
Freitag & Benz 2001), or the direct numerical solution 
of the discretized FP equation on a mesh. In this pa- 
per we focus on the latter approach. One-dimensional 
(ID) FP models assuming a spherical symmetry and 
isotropic velocity dispersion have been extensively ex- 
ploited during the last several decades, and they have 
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successfully elucidated the full dynamical history of 
star clusters (Cohn 1980; Drukier et al. 1992). Two- 
dimensional (2D) anisotropic models are a generaliza- 
tion of the ID model, in which the anisotropy of the 
velocity dispersion between the radial and the tan- 
gential directions is taken into account; however, the 
isotropized distribution function continues to be used 
for determining the diffusion coefficients (e.g., Taka- 
hashi 1995, 1996, 1997). 

The statistical methods have to make several sim- 
plifying approximations and have some limitations. 
More realistic simulations can be performed by di- 
rectly integrating the complete equations of motion 
of all stars. However, the TV-body integration still re- 
quires a huge amount of computing power. An im- 
provement in computing facilities, particularly the ad- 
vent of special purpose hardware such as GRAPE 
machines (Makino & Taiji 1998; Makino et al. 2003; 
Fukushige, Makino & Kawai 2005), made it possible 
to perform a high accuracy A r -body simulations with 
one million body (Makino & Funato 2005; Berczik et 
al. 2006; Iwasawa et al. 2006; Harfst et al. 2006). The 
comparisons between the results obtained from ID 
or 2D FP models and direct A r -body models gener- 
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ally show good agreement (e.g., Spurzem & Aarseth 
1996). Comparative studies that use the currently 
available TV-body solver for studying the dynamical 
evolution of collisional stellar systems are very impor- 
tant for checking the validity and limitations of the 
statistical methods. There have been several previous 
researches involving comparative studies (Giersz & 
Heggie 1994a,b, 1997; Giersz & Spurzem 1994; Khal- 
isi, Amaro-Seoane & Spurzem 2005; Freitag, Rasio 
& Baumgardt 2006) and the comparisons show that 
for isolated non-rotating star clusters the results of 
FP simulations are generally in good agreement with 
those of TV-body simulations. However, when a tidal 
boundary is included, a discrepancy between the N- 
body and FP models arises; this discrepancy becomes 
sensitive to TV because the relaxation and crossing 
time scales are related to different dynamical processes 
(Takahashi & Portegies Zwart 1998, 2000; Baumgardt 
2001). For example, the treatment of the tidal bound- 
ary has to be performed carefully in FP models since 
the mass evaporation process involves both orbital dy- 
namics and two-body relaxation. The TV-body model 
has to imitate the FP technique to remove stars im- 
mediately if they acquire an energy higher than the 
tidal energy; with these precautions, good agreement 
can be obtained (Spurzem et al. 2005) . 

Another extension of the ID FP model was car- 
ried out to include the effects of rotation. The first 
numerical simulation of FP models for rotating clus- 
ters was pioneered by Goodman (1983), but neither 
the results nor the code was published. A more de- 
tailed and extended work by Einsel & Spurzem (1999, 
henceforth referred as Paper I) , who developed a new 
2D FP code named "FOPAX" for this study from 
scratch, revealed that rotating clusters collapse faster 
than non-rotating ones. A post core collapse study of 
rotating star clusters, improving "FOPAX" by several 
features, including an energy source due to formation 
and hardening of three-body binaries, was done later 
by Kim et al. (2002, henceforth referred as Paper II). 
Papers I and II studied the dynamical evolution of 
rotating stellar systems using the orbit-averaged 2D 
FP equations only for equal-mass systems. In both 
cases and in the study of Goodman (1983), it was as- 
sumed that the distribution function / depends only 
on the energy and z- component of the angular mo- 
mentum (J z ). Kim et al. (2004, henceforce referred as 
Paper III) extended the method to multi-mass sys- 
tems, and exhibit interesting results concerning the 
segregation of mass and angular velocity with heavy 
stars in the cluster centres. Fiestas, Kim, & Spurzem 
(2006) modelled individual rotating globular clusters 
using the model, and Fiestas & Spurzem (2007) in- 
cluded a star-accreting central black hole with a loss 
cone. In general, according to the strong Jeans theo- 
rem there are three integrals of motion with regard to 
the orbital motion of stars in axisymmetric potentials 
(Binney & Tremaine 1987). Two of the integrals of 
motion are known as E and J z . However, the nature 



of the third integral of motion is not known clearly. 
In addition to the standard approximation required 
for FP models, the 2D FP models presented here and 
in the Papers I to III ignore any effect of the third 
integral. For stellar systems not too much flattened, 
it might be possible to construct 3D FP models fol- 
lowing the three integral models of Lupton & Gunn 
(1987), in which the square of total angular momen- 
tum (J 2 ) was adopted as an approximate third inte- 
gral. However, it would be an extremely difficult task 
numerically and physically (owing to the diffusion co- 
efficients!); even then it would still be approximate 
since we do not know the third integral analytically 
(compare also Binney & Tremaine (1987) regarding 
this issue). Hence, a thorough comparison with direct 
TV-body models appears to be the best method to ac- 
quire knowledge on the validity of all approximations 
made in our 2D FP models. 

There have been some preliminary studies of ro- 
tating star clusters involving comparisons between N- 
body and FP methods (Boily 2000; Boily & Spurzem 
2000; Ardi, Mineshige & Spurzem 2006; Ernst et al. 
2007). Although these authors found good agreements 
between these two methods, the number of cases stud- 
ied is rather limited. This can be understood in two 
ways. First, for a smaller value of TV (say up to a few 
10 4 ) a large number of statistically independent sim- 
ulations are needed, and only the ensemble average 
can be compared with the approximate FP models. 
On the other hand, different physical models, such as 
isolated or tidally limited models, different degree of 
central concentrations of stars (i.e., different central 
potential), need to be studied. In this paper, we have 
carried out a series of numerical TV-body simulations 
of rotating stellar systems, which are directly compa- 
rable with our 2D FP models in Papers I and II. By 
comparing the results with those obtained from FP 
models, we can investigate the validity of the assump- 
tions made in rotating FP models. 

This paper is organized as follows. In the next 
section, we briefly describe initial JV-body models 
and compare them with FP models. In section 3, we 
present the numerical results obtained from TV-body 
simulations and their comparisons with those of FP 
models. The summary and discussions are given in 
the last section. 



2 THE MODELS 

2.1 Numerical methods 

The Fokker-Planck (FP) code FOPAX, which takes 
into account the effect of rotation, was developed in 
Papers I, II, and III in order to study the secular evo- 
lution of star clusters having initial rotation. 

For performing direct TV-body simulations to 
study the dynamical evolution of rotating stellar sys- 
tems, we have used the currently available high ac- 
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Table 1. Initial parameters for the equal-mass models 



Wo 


u>0 


rt/r c 


rh/n 


Ttot/\W\ a 


N 




0.0 


18.0 


0.15 


0.000 


10240 


6 


0.3 


14.5 


0.18 


0.035 


10240 




0.6 


9.6 


0.24 


0.101 


10240 



a Ttot/\W\ : ratio of rotational energy to potential energy 

curacy, collisional iV-body code NBODY6 (Aarseth 
1999). The NBODY6 code uses the fourth-order Her- 
mite scheme with hierarchical block time steps (HTS) 
and the Ahmad-Cohen neighbor scheme for particle 
integration. Close encounters between stars and per- 
sistent binaries formed by three-body interactions are 
solved for their internal motion by using two-body reg- 
ularization methods (Kustaanheimo & Stiefel 1965) 
and chain regularization (Mikkola & Aarseth 1990, 
1993, 1996, 1998). Although NBODY6 is capable of 
dealing with many more astrophysical components 
such as the existence of primordial binaries and stellar 
mass-loss due to stellar evolution, we have considered 
only the treatment of close encounters between stars 
in the present study since we are mainly interested in 
the role of the initial angular momentum on cluster 
dynamics and in comparisons with FP simulations. 

2.2 Initial models and boundary condition 

In Paper II, the initial 2D FP models are generated 
according to Lupton & Gunn (1987). Our initial N- 
body models that follow rotating King models with a 
central concentration of Wo = 6 have the same con- 
ditions as those of 2D FP models in Paper II. We 
have constructed the initial models for the iV-body 
simulations from the 2D FP models in Paper II by 
random number generation. Three different initial ro- 
tations (u>o — 0.0, 0.3, and 0.6) are considered in the 
present work. In Table 1, we list some information on 
the initial models used for the present simulations. 

Since most globular clusters are bound to their 
host galaxy, they are tidally limited and stars escape 
from them through a tidal boundary. There are many 
previous studies that considered the effects of the tidal 
field on cluster evolution and comprehended it as an 
important component of the evolution (e.g., Lee & Os- 
triker 1986; Lee & Goodman 1995; Takahashi & Porte- 
gies Zwart 1998; Takahashi & Lee 2000; Yim & Lee 
2002; Lee et al. 2006; Spurzem et al. 2005). Among the 
few different implementations of modeling the tidal ef- 
fects, we adopted the instantaneous removal of stars 
whose total energy exceeded tidal energy of the clus- 
ter. This approximation is known to be inconsistent 
with the realistic iV-body treatment for small N mod- 
els, but the inconsistency decreases for a large value of 
N. We considered the equal-mass models, which are 



tidally bound to their host galaxy, in order to compare 
the obtained results with the FP results in Paper II. 
We have modified the original NBODY6 code to im- 
itate the tidal environment of the clusters modeled 
with the 2D FP equation in Paper II; this implies that 
we promptly remove stars whose energies exceed the 
tidal energy (e > end, see eq. 1), as considered in 2D 
FP models. In order to maintain the density within the 
tidal radius constant, the tidal radius decreases with 
time when there is a loss of mass through the bound- 
ary. Subsequently, we have adjusted the tidal energy 
at every regular time step according to the following 
expression, 

GM(t) 

eud (t) = (1) 

rud 

where rud and M(t) are the tidal radius and the 
total mass of the cluster within the tidal radius, re- 
spectively. We have used the initial tidal radius ob- 
tained from the FP models to compute the mean den- 
sity that is kept constant. 

The number of stars (N) in a cluster is one of the 
important parameters for the dynamical evolution of 
the cluster. While the computational burden (except 
for the core-collapse phase) does not significantly de- 
pend on N in statistical FP method, the number of 
stars is very important in the iV-body simulation as 
the computation time becomes nearly proportional to 
TV 3 . We use N = 10240 for the present equal-mass 
models only because the number should be close to 
that used in Paper II (N = 5000). In testing the va- 
lidity of our FP models, it does not matter that the 
actual number of stars in globular clusters is signifi- 
cantly larger. The choice of the number of stars deter- 
mines the relative strength of the three-body interac- 
tions, which initiate the post-collapse phase (see e.g., 
Spurzem & Aarseth 1996). 

The realization of the rotating King model for N- 
body simulations is shown in Figs. 1 and 2. We have 
shown the radial profiles of density for both the N- 
body and FP models (Fig. 1) and the distribution of 
the radial and tangential velocities of the stars in the 
iV-body realization of the rotating King models (Fig. 
2). Three iV-body models having different degrees of 
rotation are compared with the 2D FP models. Each 
density profile is obtained from the mean of 10 dif- 
ferent initial models generated by different random 
seed numbers. The open circles represent the density 
profile adopted in the FP models. Each rtid,o is the 
initial tidal radius derived from the FP model. Excel- 
lent agreement is observed between the radial profile 
beyond the core radius (r c ) in the iV-body realization 
and that in the FP model. However, within r c , the 
density of the iV-body model is slightly lower than 
that of FP model. This may be due to the fact that the 
number of stars inside the core is rather small. How- 
ever, due to inevitable random fluctuations, it is im- 
possible to construct initial iV-body models perfectly 
identical to the FP models. We believe that statistical 
FP models agree very well with the averaged iV-body 
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Figure 1. Radial profiles of the density for the initial N- 
body models with a central potential of Wo = 6 and loq = 
0.0, 0.3, and 0.6. The upper-left panel shows the FOPAX 
results. For comparison, the density profile of the initial 
FP model with (W ,w ) = (6, 0.6), (6, 0.3), and (6,0.0) is 
shown (open circles). The JV-body realization of the initial 
model shows good agreement with that of the FP model, 
except for the central region. 



Figure 2. Tangential and radial velocity distributions 
of the initial TV-body models with a central potential of 
Wo = 6 and initial rotation of u>o = 0.0 (left panels) and 
0.6 (right panels). The X-axis represents the sky-projected 
equatorial distance measured in units of initial core radius. 
The mean radial velocity distribution along the projected 
equator is shown by filled circles with lu errors. The central 
solid body rotation and the subsequent highly differential 
rotation are typical of rotating clusters (see Fig. 10 of Pa- 
per II for comparison). 



3 RESULTS 



models and increasing the number of stars will im- 
prove the degree of agreement. Fig. 1 also shows that 
increasing the rotation decreases the concentration of 
the stellar system (smaller ratio between the tidal and 
the core radii). 

The position-velocity distributions that are sky- 
projected are displayed in Fig. 2 for the iV-body mod- 
els of non-rotating (a>o = 0.0, left panels) and highly 
rotating (luo = 0.6, right panels) clusters with the cen- 
tral potential of Wo = 6. To have maximum effect of 
intial rotation in sky projectied distribution we project 
the model clusters on sky in such a way that the ro- 
tating axis is a perpendicular axis in sky-projected 
plane. The distance to the rotation axis is measured 
in units of the initial core radius. We also show the 
mean profiles of the radial velocities with la errors 
(filled circles). The profile of the radial velocity for 
the non-rotating model shows increasing velocity dis- 
persion toward the center of the cluster. A typical ve- 
locity structure of the rotating model is shown clearly 
in the distribution of the radial velocities of the rotat- 
ing model: rigid-body rotation inside the core radius 
and highly differential rotation subsequently. 



3.1 Core collapse, central density, and 
central velocity dispersion 

We start the discussions by presenting the results for 
the evolution of the central density and central veloc- 
ity dispersion (Fig. 3). The central density increases 
with time due to the two-body relaxation for an equal- 
mass system, where the time is measured in units of 
initial half-mass relaxation time {r r h,o)- The expres- 
sion for T r h.o for the equal-mass system is given by the 
following formula (Spitzer & Hart 1971): 



Trhfi = 0.138 



G l/2 m l/2 lnA ' 



(2) 



where N, r hy0 , G, m, and In A = In^N) denote the 
total number of stars, initial half-mass radius, grav- 
itational constant, mean mass of stars and Coulomb 
logarithm, respectively. It has been shown by Giersz & 
Spurzem (1994) and Giersz & Heggie (1994a,b) that 
the best agreement between the direct A-body cal- 
culations and orbit-averaged FP equation is achieved 
when the coefficient 7 in the Coulomb logarithm has 
a value of 0.11. Therefore we use this value in the 
present work as it is also used in Paper II. Since the 
half-mass radius varies with the rotation parameter 
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Figure 3. Evolution of central density and central velocity 
dispersion for each models. FOPAX models arc shown with 
smooth lines. 
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Figure 4. Time evolutions of the Lagrangian radii contain- 
ing 1%, 5%, 10%, 20%, 50%, and 75% of the initial cluster 
mass. The initial parameters that characterize the cluster 
model are written on the lower-left corner of each panel. 
The FOPAX calculations are shown with dashed lines on 
each panel. 

loo, the values of T r h,o could also depend on coo, even 
for models with the same Wo. 

The time of core collapse (t cc ), the time for the 
complete disruption of the cluster (tdis), and the clus- 
ter mass at the time of core collapse in units of initial 
cluster mass (M cc ) are listed in Table 2. We can easily 



notice that the rotating models evolve faster than the 
non-rotating ones, in both the Af-body and FP mod- 
els. The faster the cluster rotates the shorter the time 
taken for the core collapse. As discussed in detail in 
Papers I & II, the acceleration is caused by the com- 
bination of gravothermal and gravogyro instabilities. 

In Fig. 3, we show the evolution of the central 
density (p c ) and central velocity dispersion (er c ) of the 
cluster. The results obtained from the FP and A'-body 
models are displayed as smooth lines and solid lines 
with a large fluctuation, respectively. Although the 
core collapse occurs at slightly different times (one 
can estimate t cc more accurately in Fig. 4), there 
is good agreement between the FP and A'-body re- 
sults. The A'-body model produces rather noisy data 
since there is a significant statistical fluctuation in 
the physical parameters. However, we can still per- 
form some quantitative comparisons between the FP 
and the A-body approaches. In the early evolutionary 
stage, the central density derived from the A'-body 
is less than the value obtained by FP for the model 
with ujq = 0.6. This may reflect the difficulty in de- 
termination of the central density for Af-body models, 
especially for rapidly rotating clusters that are signifi- 
cantly flattened. Therefore, the most rapidly rotating 
model has the largest discrepancy with regard to the 
central density between the N-body and the FP mod- 
els. 

From Table 2, we notice that the times for core 
collapse in the FP simulations are generally greater 
than those in the N-body simulations by 5-15%. Ap- 
parently, different approaches should yield some dif- 
ferent results, and we regard a small difference of 5- 
15% in t cc as being insignificant. These differences 
would decrease for large N models as the assumptions 
made in the FP method become more appropriate. In 
fact, we have observed that t cc decreases for model 
with large A" value when we compared the simula- 
tions performed for A-body models with N = 5000 
and N = 10240. 

After the core collapse, the evolutions of p c de- 
rived from the A'-body and FP models are also some- 
what different from each other. Toward the end of 
the evolution, the difference becomes quite significant. 
The disruption time of the A'-body models is slightly 
larger than that of the FP models. This is due to the 
fact that the escape rate varies with the number of 
stars. We perform more detailed comparisons on the 
evaporation of stars in section 3.3. It is easier to deter- 
mine the exact core collapse time based on the evo- 
lution of the central velocity dispersion because the 
fluctuation amplitude is smaller than that of the cen- 
tral density. The times of core collapse listed in Table 
2 are determined by inspecting the behavior of a c . 

We show the evolution of the Lagrangian radii 
of the equal-mass models in Fig. 4. The results from 
the FP simulation are displayed by the dashed lines. 
Each line represents the radii where the cluster con- 
tains 1%, 5%, 10%, 20%, 50%, and 75% of the initial 
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Table 2. Time scales of tidally bound models. 



Model 


Wo 




tcc[r r h,o] 


tdis[Trhfi] 


M cc 






0.0 


10.1 


24.30 


0.63 


JV-body 


6 


0.3 


8.7 


17.75 


0.55 






0.6 


6.9 


11.55 


0.37 






0.0 


11.73 


22.61 


0.59 


FOPAX 


6 


0.3 


10.31 


16.96 


0.48 






0.6 


7.27 


10.08 


0.33 
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Figure 5. The evolution of a c as a function of p c . a c fol- 
lows power laws during pre- and post-collapse states. The 
evolution is nearly independent of luq . 



mass of the cluster. It is not straightforward to de- 
termine the Lagrangian radii in a flattened system. 
In Paper I, the Lagrangian radii were evaluated along 
the specific zenith angle where the effects of flattening 
on the mass shells are expected to be less important; 
the same zenith angle is used in Fig. 4. However, as 
our models are nearly spherical, we determine the La- 
grangian radii on the assumption that the system is 
spherically symmetric for the JV-body models. As seen 
in Fig. 4, differences between the FP and the A^-body 
models are very small. Analyzing in more depth, we 
notice that after the core collapse, the inner part of the 
cluster expands more rapidly in the FP model than in 
the A-body model, although the difference is rather 
small. 

The relationship between the central density and 
the central velocity dispersion is shown in Fig. 5. The 
upper-left panel (Fig. 5a) shows the relation between 
a c and p c which is obtained from all the three FP 
models with different initial rotations. The other three 
panels (Figs. 5b, 5c, and 5d) show the results of the N- 



body models. From Papers I and II, we know that this 
relationship is not affected considerably by the initial 
rotation, as shown in Fig. 5a. Again, we find good 
agreement between the FP and the A-body models. 
The large amount of scatters during the post core col- 
lapse phase, as shown in Fig. 5 is mainly due to the 
large fluctuation in the central density. The power-law 
behavior of a c on p c during the pre-collapse phase is 
a consequence of the self-similarity of the collapsing 
core and it is well known that <7 C oc p c during this 
stage (Cohn 1980). During the post-collapse phase, 
we can derive the relationship between a c and p c us- 
ing an energy balance argument and the assumption 
of self-similar evolution (see section 3.4). It follows 
that a c oc with /3 = 0.25, which is in good agree- 
ment with P — 0.23 derived from the FP model in 
Paper II. The A^-body results appear to follow similar 
power-law behavior, although the power-law index (5 
is difficult to determine because of the large scatter. 

From the A^-body calculations, we have confirmed 
the earlier finding of a significant acceleration in clus- 
ter evolution due to rotation, which was obtained from 
the FP calculations in Papers I and II. We also find 
that both the A-body and FP models give similar re- 
sults, although there are small differences with regard 
to the time of core collapse, and the disruption times. 



3.2 Evolutions of anisotropy and angular 
momentum 

In axially symmetric systems, the natural decomposi- 
tion of velocity vectors is to use the cylindrical coordi- 
nate which has its origin at the center of mass of the 
cluster. We investigate the evolution of the velocity 
dispersions (an, a^, and a z ) and show the evolution 
of these quantities in Fig. 6, where (R, 4>, z) represents 
the conventional axis of the cylindrical coordinate sys- 
tem. 

For initially rotating models, these three velocity 
dispersions have different values. In the right panels of 
Fig. 6, the ratio of to oo is initially greater than 1 
for the rotating models and approaches the isotropic 
value of 1, where <to represents the average ID velocity 
dispersion defined by cto = (cfj + o\ + vl) 1 ^ 2 /V3. 

The angular momentum in rotating stellar sys- 
tems is transferred outward through two-body relax- 
ation, and is also lost due to escaping stars. The stars 
gaining a large angular momentum migrate to the 
outer parts of the system, while those losing angu- 
lar momentum drift toward the central parts. As the 
stars with a high angular momentum move outward 
and finally escape from the cluster, the total angular 
momentum of the system decreases with time. 

In Fig. 7, we display the evolutions of the z- 
component of the angular momentum per unit mass 
(J z ) of an entire cluster. The results from the A"-body 
models are shown by solid lines, while those from the 
FP models are shown by dashed lines. To indicate the 
degree of error in determining J z , we also show the 
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Figure 6. Time evolutions of 072, a^, and a z (left panels) 
and a^ /<To (right panels) for each models. At the time of 
core collapse, an , a§ , and az are almost equal and a<p / a® 
is Ri 1, where ao is the ID velocity dispersion. 
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Figure 7. Evolution of the z-componcnt of angular mo- 
mentum per unit mass (J z ) (Fig. 7a) and the loss rate of 
J z (Fig. 7b). The results for the JV-body model and FP 
model arc shown with solid lines and dashed lines, respec- 
tively. The J z evolution for the non-rotating model indi- 
cates the corresponding error in determining the J z . It is 
clearly shown that J z for rotating models decreases with 
time due to the escape of stars possessing an angular mo- 
mentum and that the time evolutions of J z in the N-hody 
and FP models agrees quite well. The loss rate of J z for the 
./V-body models (solid lines) is smoothed for easier compar- 
ison. 



Figure 8. Distribution of the at 4 selected evolutionary 
epochs in cylindrical coordinate system for a model with 
{Wq,u>q) = (6,0.6). When there is no rotation, it should 
show a symmetry with respect to = 0. A asymmetry of 
V0 disappears with increasing time due to the loss of the 
angular momentum. The filled circles are the average 
values and the error bar corresponds to 1 a dispersion. 

time evolution of J z for the non-rotating model. We 
first note that there is good agreement between the N- 
body and FP results. The effect of neglecting the third 
integral in the FP models is minimized when dealing 
with global cluster properties (e.g., total cluster mass 
and total angular momentum), while this effect is se- 
vere for local properties (e.g., central density). It is 
clearly shown that J z monotonically decreases with 
time. The loss rate of angular momentum is large in 
the early phases and decreases as the cluster evolves 
(see Fig. 7b). The combined effect of gravitation and 
rotation accelerates the evolution of the cluster. A sub- 
stantial loss of the initial angular momentum in an en- 
tire cluster prevents rotation from playing an impor- 
tant role in the evolution of a cluster in later phases. 
Since the cluster is losing mass at a nearly constant 
rate, the total angular momentum of the cluster de- 
creases more rapidly than J z . 

3.3 Evolutions of rotation curve 

We can construct the rotation curves by computing 
the averages of V^, which is the ^-component of the 
velocity of stars. We choose the model with loo = 0.6 
and present the distribution of V$ for all stars at some 
specific epochs in Fig. 8. The asymmetry of with 
respect to V$ = indicates the global rotation of the 
stellar system. During the pre-collapse phase, V$ be- 
comes more dispersive, particularly around the cen- 
tral region. A large spread near the center (R ~ 0) at 
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Figure 9. Radial profiles of mean rotational velocities of 
the cluster model with (Wo,u>o) = (6,0.6) at four evolu- 
tionary stages. The radii are measured in units of current 
core radius (r c ). For comparison, we show the rotational 
profiles from the FP models by dashed lines. 




Figure 10. Density distribution of stars in (R, Z) coordi- 
nate at four selected evolutionary epochs for a model with 
(Woj^o) = (6,0.6) as computed by N-body (contour map) 
and FP methods (color map). The horizontal axis is R run- 
ning from to 3 in units of the initial core radius, while 
the vertical axis is the z coordinate with the same scale. 
The ./V-body and FP methods give almost the same density 
distribution. 



t = t cc and the slight asymmetry of V$ indicate that 
the cluster is losing angular momentum; this also in- 
dicates that as the core collapse approaches, the stars 
are falling into the central region of the cluster with a 
high angular speed (but low rotation velocity). After 
tcc, V<j> still has a small amount of asymmetry. 

The radial profiles of the rotational velocity (tVot) 
for the model with ujq = 0.6 are shown in Fig. 9. We 
also display the radial profiles of v ro t obtained from 
the FP model by dashed lines. The profiles at different 
evolutionary epochs are shown with different symbols: 
filled squares for the initial model, open circles for the 
pre-collapse phase, filled circles at the time of core 
collapse, and open star marks for the post-collapse 
phase. The radial profiles in Fig. 9 are obtained by 
averaging for stars having the same value of R in 
cylindrical coordinates. 

The radial profiles of v ro t at the initial epoch 
shows good agreement between the iV-body and FP 
models except for a very central region (R < 0.5r c ). 
The amount of rotation decreases with time and the 
radius where v ro t has the peak value progresses out- 
ward when the radius is measured in units of current 
core radius. This is a consequence of the self-similar 
evolution during the pre-collapse phase (Papers I & 
II). The radial profiles of v ro t for the TV-body model 
during pre-collapse phaes agree well with those for the 
FP model. However, at later times, the two methods 
yield somewhat different rotation curves. For example, 
at the time of core collapse, the TV-body model pre- 
dicts a larger rotation speed for the stars in the cen- 
tral and outer parts as compared to the FP method. 
The rotation seems to continue until a very late epoch 
(t = 10t r h,o) in the TV-body model, while there is vir- 
tually no rotation left in the FP model. 

In Fig. 10, we have shown the distribution of the 
stellar density in the form of a color map for the 
FP results and in the form of isodensity contours for 
the iV-body results. The four panels represent differ- 
ent epochs for the cluster with largest initial rotation 
speed (uio — 0.6). We display the cluster shapes in 
(R, z) coordinates because we assume that the initial 
model has the axis of rotation along z in cylindrical 
coordinates. To compare the cluster shapes at differ- 
ent epochs, we measure the size in units of initial core 
radius (r Cj o). First, we notice that the shapes com- 
puted by the iV-body and FP models are almost iden- 
tical. As the rotation becomes negligible, the stellar 
system becomes more spherical. The effects of rotation 
on the flattened shape of a cluster are observed during 
the time of core collapse and the post-collapse phase 
(Figs. 10c & lOd) ; the shape of the cluster around the 
central region is almost spherical. 



3.4 Tidal boundary 



If the cluster rotates around the host galaxy on a cir- 
cular orbit, the tidal field experienced by the cluster 
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Figure 11. The evolution of the total mass of the cluster. 
The discrepancy between FOPAX and NBODY6 is mainly 
due to the number of particles after the core collapse. 

does not change with time. In a steady tidal field, the 
tidal radius is expressed as follows: 

/ M \ 1/3 

r ^\m-o) Rg > (3) 

where M is the total mass of the cluster within the 
tidal boundary; Rg, the distance of the cluster from 
the galactic center; and Mg, the galactic mass within 
Rg- The above equation ensures that the mean den- 
sity within the tidal radius is a constant. As the cluster 
loses the stars beyond the tidal boundary, the tidal ra- 
dius has been adjusted to maintain a constant mean 
density. 

We depict the evolution of the total mass of the 
cluster in Fig. 11. During the pre-collapse phase, the 
evolution of the total mass, which is computed from 
the A-body and FP models, are similar but shows 
small deviation each other. However, after the core 
collapse, the two methods result in somewhat differ- 
ent bevaiour. This descrepancy in post-collapse phase 
is partly due to the accumulation of small difference 
during the pre-collapse phase. The escape rate of the 
A-body simulation is lower than that of the FP simu- 
lation. This may be caused due to the decrease in the 
number of stars in the cluster. We have assumed the 
instantaneous escape of stars whose energy exceeds 
the tidal energy. However, in the FP simulations, the 
time step is determined by the fraction of the relax- 
ation time while the A-body uses a time step pro- 
portional to the crossing time scale. As the cluster 
loses its mass, the ratio between these two time steps 
changes. In other words, the crossing time scale and 
the relaxation time scale have the following relation- 
ship (Binney & Tremaine 1987): 



Figure 12. Time evolution of r c , Th, and rud measured 
in units of initial half-mass radius (Yj, q). The most rapidly 
rotating cluster has the largest r c and r^, and the smallest 
rtid- 



irelax ~ ^ Across- (4) 

Therefore, t cr oss/t re i a x increases as A decreases. This 
means that the A-body simulation removes stars less 
effectively than the FP simulation for small value of 
A; this is why we observe a long tail in the A-body re- 
sults in Fig. 11. In A-body simulations, after the core 
collapse, the number of stars remaining in the cluster 
is of the order of 10 3 . This is not a sufficient number 
to obtain results comparable with the FP method. 
Hence, the escape rate for A-body models is lesser 
than that for the FP models during the late stages of 
the evolution and A-body models survive longer than 
FP models. We can also observe this effect in Fig. 3. 
In this figure, central density features are observed to 
be inconsistent between the A-body and the FP mod- 
els toward the end of the evolution. With more stars, 
this gap would become narrower. Since the number of 
stars in a real globular cluster is considerably larger 
than that used in the present A-body simulations, the 
difference between A-body and FP models would de- 
crease for a realistic number of stars (see Takahashi & 
Portegies-Zwart 1998). 

3.5 Core, half-mass, and tidal radii 

We now investigate the evolution of the core, half- 
mass, and tidal radii of the star cluster. The behavior 
of r t id and r c in units of initial half-mass radius (r%o) 
is shown in Fig. 12. As the rapidly rotating initial 
model has a half-mass radius smaller than the slowly 
rotating model, the initial value of r c /rh for the most 
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pj/ 6 'M 1 / 3 . On the other hand, according to Goodman 
(1987), the energy balance argument predicts that 
pel ph oc M 4 / 3 . Therefore, o c ~ p l J A p^ 1 ^ 12 . II we use 
the tidal boundary condition, M/r 3 id = constant and 
Ph ~ t- , we can obtain the following relation: 

In <r c ~ — In p c — — In . (5) 

4 Fc 4 r h K ' 

The variation in during the post-collapse phase 
is very small (by a factor of few) as compared with 
that in p c (by a few orders of magnitude), except near 
the disruption time. Therefore, we can approximately 
write following relation: 

lna c ~ i lnp c . (6) 

If we express a c oc pf , j3 — 0.25, then. This value is 
close to 0.23 that was achieved in Paper II. 



Figure 13. Time evolution of r c /rh and rtid/rh- Self- 
similarity feature after t cc is shown. In particular, the evo- 
lution of r t id/rh shows good agreement with that of FP. 

rapidly rotating cluster is larger than that of the other 
two models, as shown in Fig. 13 (Paper II). The evo- 
lutions of r c and rud in units of rh at the same evolu- 
tionary stage as those in Fig. 12 are shown in Fig. 13. 
We find that there are significant differences in r c /rh 
between the FP and the iV-body models, while the 
difference is not so apparent in r tid /r h . This reflects 
the difficulties in determining r c for the iV-body mod- 
els with a relatively small value of N rather than any 
systematic differences in different approaches. Both r c 
and rh decrease with time, although there is a differ- 
ence in the decreasing rate that depends on the initial 
degree of rotation. After the core collapse, both r c 
and rh increase for some time at almost the same rate 
due to self-similar expansion. Subsequently, the tidal 
boundary shrinks rapidly after core bounce. There- 
fore, the half-mass radius begins to decrease again. 
However, the shrinking of r U d does not affect the cen- 
tral region and r c continues to increase. At the end 
of the evolution r c shows a steep increase; this signals 
the complete disruption of the cluster. 

After the core collapse, r c /rh and r ti d/rh show 
nearly the same behavior. The value of r c /rh is al- 
most constant and the evolution of rud/rh for different 
initial models is similar to each other and indepen- 
dent of the rotation parameter ujo- We already have 
assumed the self-similar evolution of the inner part 
of the stellar system in order to explain the relation 
between p c and a c after t cc (Fig. 4). We can express 
°"c ~ ^r- and °~'h ~ ~^ tL , where M c and Mh are the 
masses within r c and rh, respectively. Since the inner 
parts of the cluster are nearly isothermal, we obtain 
o"h/cr c = constant. With the self-similarity assump- 
tion (r c /rh = constant), we can rewrite a c as o c ~ 



4 SUMMARY AND DISCUSSION 

We have performed numerical simulations for the evo- 
lution of initially rotating star clusters with equal 
mass using NBODY6 and have compared the results 
with those computed by the direct integration of the 
Fokker-Planck equation. We have considered clusters 
with N = 10240. For critical comparisons between 
A"-body and FP models, we constructed the initial N- 
body models using the initial 2D distribution function 
used for FP models. 

We observed the acceleration of the core collapse, 
as reported in Papers I & II. The degree of accel- 
eration obtained from the present iV-body models is 
slightly different from that obtained from the FP mod- 
els; however, the small difference in the core-collapse 
time between the iV-body and statistical methods (FP 
model, gaseous model, etc.) has also been observed 
earlier (Spurzem & Aarseth 1996). The entire evolu- 
tionary trend of the central density agrees with that 
of the FP models. 

The ^-component of the specific angular momen- 
tum (J z ) is observed to monotonically decrease with 
time for the clusters with initial rotation. The global 
evolutionary trend of J z between the iV-body and the 
FP models shows excellent agreement. The loss rate 
of J z decreases as the cluster evolves. Therefore, we 
conclude that during the early stages the existence of 
initial rotation significantly affects the entire cluster 
evolution. 

In FP simulations, the cluster evolution will be 
independent of the third integral to the end of time. 
On the other hand, in the iV-body simulations, the 
third integral effect may appear during the evolution. 
In addition, there is a limit on the number of stars and 
the random fluctuations of the iV-body models in this 
study and these limits also causes differences with the 
FP method. Therefore, we need to perform more N- 



© 0000 RAS, MNRAS 000, 000-000 



Comparative study between N-body and Fokker-Planck simulations for rotating star clusters: I. Equal-mass s 



body simulations with a larger number of stars or with 
various models by using different random number. 
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